home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / 3d-tools / dog2.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-06-06  |  16.0 KB  |  608 lines
  1. /*
  2. * dog.c 2.0 -    filter an image by applying difference of Gaussians mask.
  3. *        The input is in byte or float format, and the output is in
  4. *        floating point.    or in integer format.            */
  5. char    usage[]="\n\
  6. * usage:    dog [-A esigma [masksize [ratio ]]] [-b #] [-n #]    \n\
  7. *        [-p #] [-i [-c]] [-m] [-g #] [-w] [< input_seq] [> out_seq]\n\
  8. *\n\
  9. *    esigma [1.0] is the standard deviation of the \"excitatory\" Gaussian.\n\
  10. *        A real positive number.    \n\
  11. *    masksize [7] is the size of the mask (an integer).        \n\
  12. *    ratio [1.6] is the ratio b/w the s.d.'s of the inhibitory and    \n\
  13. *        excitatory Gaussians.    \n\
  14. *    -p followed by a positive integer specifies the precision.    \n\
  15. *        Defaults to 1.    \n\
  16. *    -b begin from #th frame. default is 0.    \n\
  17. *    -n number of frames will be processed, default is total.    \n\
  18. *    -i implies output in IFMT_LONG format.    \n\
  19. *    -c if -i is specified, causes checking of input to be in the    \n\
  20. *     range [-1024 to 1024].    \n\
  21. *    -m output the Gaussian(s) only, do not convolve.        \n\
  22. *    -g display size = filter size + #    \n\
  23. *        otherwise, Gaussian filter display by using 2**n square.\n\
  24. *    -w output without header.    \n\
  25. *\n\
  26. *    If input file is not given, the filtering is done on an impulse    \n\
  27. *        response in a 7 x 7 frame.    \n";
  28. /*
  29. * compile:    cc -O -o DestPath/dog dog.c -lccs -lhips -lm
  30. *  note:    The efficiency of convolution depends heavily on optimization
  31. *    by    the compiler, hence option -O, and also need to define Faster.
  32. *
  33. * Author:    Yoav Cohen - 12/12/82
  34. * modified by Jin, Guojun for convolution algrithm correction - 10/25/90
  35. */
  36.  
  37.  
  38. #include "header.def"
  39. #include "imagedef.h"
  40. #include <math.h>
  41.  
  42. U_IMAGE    uimg;
  43.  
  44. #define    rows    uimg.height
  45. #define    cols    uimg.width
  46. #define    frm    uimg.frames
  47. #define    k    num_offp
  48. #define    bbuf    uimg.src
  49. #define    fname    uimg.name
  50.  
  51. #define    Float    float
  52. #define    Range    1024
  53. #define    Faster
  54.  
  55. #ifndef    MType
  56. #define    MType    int
  57. #endif
  58.  
  59.  
  60. main(argc, argv)
  61. int    argc;
  62. char**    argv;
  63. {
  64. MType    *int_inbp, *int_obp, *int_conv1,
  65.     *int_emask, *int_imask, fsize;
  66. int    fr, synthetic=0,
  67.     woh=0, bgn_frm=1, num_offp=0,
  68.     isint=0, checking=0, gaussonly=0,
  69.     precision=1,
  70.     masksize=7,
  71.     i, edglen,
  72.     ovf, unf;
  73.  
  74. #define    form_size    gaussonly
  75. #define    c    edglen
  76.     /*    using same variable at different parse    */
  77.  
  78. char    *bp;    /* byte buffer & pointer    */
  79. Float    *inbuf, *obuf,    /* buffers    */
  80.     *conv1,        /* temp converting buffer    */
  81.     *obp, *inbp,    /* in & out put buffer pointers    */
  82.     gaussize=0,
  83.     *emask, *imask;
  84. double    esigma=1., ratio=1.6,    /* default value    */
  85.     gauss_mask();
  86.  
  87.  
  88. format_init(&uimg, IMAGE_INIT_TYPE, HIPS, -1, *argv, "S20-1");
  89.  
  90. for (i = 1; i < argc; i++)
  91.   if(*argv[i] == '-')
  92.     {    c = 1;
  93. args:    switch(argv[i][c++])
  94.     {
  95.     case 'p':
  96.         if (!argv[i][c]) {    i++; c=0;    }
  97.         precision = atoi(argv[i]+c);    break;
  98.     case 'A':
  99.         if (!argv[i][c]) {    i++; c=0;    }
  100.         esigma    = atof(argv[i++]+c);
  101.         if (i >= argc)    break;
  102.         if (*argv[i] == '-')    goto    args;
  103.         frm = masksize= atoi(argv[i++]);
  104.         if (i >= argc)    break;
  105.         if (*argv[i] == '-')    goto    args;
  106.         ratio    = atof(argv[i]);    break;
  107.     case 'i':
  108.         isint++;    break;
  109.     case 'c':
  110.         checking++;    break;
  111.     case 'b':
  112.         if (!argv[i][c]) {    i++; c=0;    }
  113.         bgn_frm = atoi(argv[i]+c);    break;
  114.     case 'n':
  115.         if (!argv[i][c]) {    i++; c=0;    }
  116.         num_offp = atoi(argv[i]+c);    break;
  117. #ifdef    _DEBUG_
  118.     case 'd':
  119.         debug++;    break;
  120. #endif    _DEBUG_
  121.     case 'm':
  122.         gaussonly++;
  123.         msg("%s no convolution; output: the Gaussian(s)\n", Progname);
  124.         break;
  125.     case 'g':
  126.         if (!argv[i][c]) {    i++; c=0;    }
  127.         gaussize=atof(argv[i]+c);    break;
  128.     case 'w':
  129.         woh++;    break;
  130.     default:
  131. info:        usage_n_options(usage, i, argv[i]);
  132.     }
  133.    }
  134.    else    fname = argv[i];
  135.  
  136. io_test(fileno(out_fp), goto    info);
  137.  
  138. if (esigma <= 0.0)    prgmerr('e', "sigma must be positive");
  139. if (masksize < 1)    prgmerr('m', "mask-size must be > 1");
  140. if (ratio < 0.0)    prgmerr('r', "ratio must be >= 0");
  141. if (ratio == 1)        ratio = .99995;
  142.  
  143. if (checking && !isint) {
  144.     mesg("checking mode ignored; output is not INT.\n");
  145.     checking=0;
  146.     }
  147.  
  148. io_test(fileno(in_fp), synthetic= !iset);
  149. if (fname)    {
  150.     in_fp = freopen(fname, "rb", stdin);
  151.     if (synthetic = !in_fp)    syserr("input file %s", fname);
  152.     goto    rdh;
  153. }
  154. else if (synthetic)    {    /* no input file */
  155.     edglen = 1;
  156.     if (gaussize)
  157.         edglen = masksize + gaussize;
  158.     else    while(edglen <= masksize)    edglen <<= 1;
  159.     if (!gaussonly)
  160.         fprintf(stderr,"%s: output frame size = %dX%d.\n",
  161.             Progname, edglen, edglen);
  162.     rows = cols = edglen;
  163.     uimg.pxl_out = sizeof(float);
  164.     uimg.o_form = IFMT_FLOAT;
  165.     uimg.frames = 1;
  166.     (*uimg.std_swif)(FI_INIT_NAME, &uimg, "DOG mask", 0);
  167. } else    {
  168. rdh:    (*uimg.header_handle)(HEADER_READ, &uimg, 0, 0);
  169.     if (uimg.in_form != IFMT_BYTE && uimg.in_form !=IFMT_FLOAT)
  170.         prgmerr(1, "image pixel format must be bytes or floats");
  171.     }
  172.  
  173. fsize = rows*cols;    /* size of one frame */
  174.  
  175. /* convert number of frames for process to last frame number,
  176.     it has different meaning from 3d_dog.
  177. */
  178. if (!bgn_frm || bgn_frm>frm)
  179.     bgn_frm=0;
  180. else    bgn_frm--;
  181. if (num_offp && num_offp+bgn_frm < frm)    frm = num_offp + bgn_frm;
  182.  
  183. if (!ratio)    message("%s: blurring by a Gaussian mask\n", Progname);
  184.  
  185. if (!gaussonly)    {
  186.     if (uimg.in_form==IFMT_BYTE)
  187.         bbuf = nzalloc(fsize, 1L);
  188.     inbuf = nzalloc(fsize, (MType)sizeof(Float));
  189.     conv1 = nzalloc(fsize, (MType)sizeof(Float));
  190.     obuf = nzalloc(fsize, (MType)sizeof(Float));
  191.     uimg.o_form = IFMT_FLOAT;
  192.     if (isint) {
  193.         int_inbp = (MType *)inbuf;
  194.         int_conv1 = (MType *)conv1;
  195.         int_obp    = (MType *)obuf;
  196.         int_emask = zalloc((MType)masksize, (MType)sizeof(MType));
  197.         int_imask = zalloc((MType)masksize, (MType)sizeof(MType));
  198.         uimg.o_form = IFMT_LONG;
  199.     }
  200.     uimg.pxl_out = 4;
  201.     if (!woh)
  202.         (*uimg.header_handle)(HEADER_WRITE, &uimg, argc, argv, True);
  203.     else    (*uimg.header_handle)(HEADER_FWRITE, &uimg, argc, argv, stderr);
  204. }
  205.  
  206. emask = zalloc((MType)masksize, (MType)sizeof(Float));
  207. imask = zalloc((MType)masksize, (MType)sizeof(Float));
  208. gauss_mask(esigma,masksize,emask,precision,gaussonly?out_fp:NULL);
  209. if (ratio)
  210.     gauss_mask(esigma*ratio,masksize,imask,precision,gaussonly?out_fp:NULL);
  211.  
  212. if (gaussonly)    return(0);    /* Done    */
  213.  
  214.     if (isint)
  215.     {  for(i=0; i < masksize; i++) {
  216.         int_emask[i] = emask[i] * Range + .5;
  217.         int_imask[i] = imask[i] * Range + .5;
  218.         }
  219.         form_size = sizeof(MType);
  220.     }
  221.     else    form_size = sizeof(Float);
  222.  
  223. if (bgn_frm)    fseek(in_fp, (bgn_frm)*fsize, 1);
  224.  
  225. for(fr = bgn_frm; fr < frm; fr++)
  226.    {    if (frm > 1)
  227.         fprintf(stderr,"%s: starting frame #%d\n", Progname, fr);
  228.     if (synthetic) {
  229.         if (edglen & 1)
  230.             synthetic = (rows*cols >> 1) + 1;
  231.         else    synthetic = (rows*cols + cols) >> 1;
  232.         if (isint) {
  233.             for(i=0; i < fsize; i++)    int_inbp[i] = 0;
  234.             int_inbp[synthetic] = 1; /* center = 1 */
  235.         }
  236.         else {
  237.             for(i=0; i < fsize; i++)    inbuf[i] = 0.0;
  238.             inbuf[synthetic] = 1.0;
  239.         }
  240.     }
  241.     else{    /* read one frame image --> memory    */
  242.         if (uimg.in_form==IFMT_BYTE)
  243.            {    (*uimg.std_swif)(FI_LOAD_FILE, &uimg, uimg.load_all=0, No);
  244.             if (isint)    /* type convert    */
  245.                for(i=0, bp=bbuf, int_inbp=(MType*)inbuf; i<fsize; i++)
  246.                 *int_inbp++ = *bp++ & 0xFF;
  247.             else for(i=0, bp=bbuf, inbp=inbuf; i<fsize; i++)
  248.                 *inbp++ = *bp++ & 0xFF;
  249.            }
  250.         else{
  251.             if (upread(inbuf, form_size, fsize, in_fp) != fsize)
  252.                 syserr("unexpected EOF");
  253.             if (isint)     /* type convert; Float -> int */
  254.                for(i=0,int_inbp=(MType*)(inbp=inbuf); i<fsize; i++)
  255.                 *int_inbp++ = *inbp++;
  256.             }
  257.         }
  258.  
  259.     if (checking)
  260.        for(i=unf=ovf=0,int_inbp=(MType*)inbuf; i<fsize; i++){
  261.         k = *int_inbp++;
  262.         unf += (k < -Range);
  263.         ovf += (k > Range);
  264.        }
  265.     /*    starting convolution    */
  266.     if (isint) {
  267.         int_vconvolve(int_emask,masksize,int_inbp,int_conv1,rows,cols);
  268.         int_hconvolve(int_emask,masksize,int_conv1,int_obp,rows,cols);
  269.         if (ratio!=0.0) {
  270.             int_vconvolve(int_imask,masksize,int_inbp,int_conv1,rows,cols);
  271.             int_hconvolve(int_imask,masksize,int_conv1,int_inbp,rows,cols);
  272.         for(i=0, int_obp=(MType*)obuf, int_inbp=(MType*)inbuf; i<fsize; i++)
  273.             *int_obp++ -= *int_inbp++;
  274.         }
  275.     }
  276.     else {
  277.         vconvolve(emask,masksize,inbuf,conv1,rows,cols);
  278.         hconvolve(emask,masksize,conv1,obuf,rows,cols);
  279.         if (ratio) {
  280.             vconvolve(imask,masksize,inbuf,conv1,rows,cols);
  281.             hconvolve(imask,masksize,conv1,inbuf,rows,cols);
  282.             for(i=0, obp=obuf, inbp=inbuf; i < fsize; i++)
  283.                 *obp++ -= *inbp++;
  284.         }
  285.     }    /* end of convolution    */
  286.     fwrite(obuf, form_size, fsize, out_fp);
  287.    }/*    end for frames    */
  288.  
  289. if (checking && (ovf || unf))
  290.     prgmerr(0, "%d overflows or %d underflows !!\n", ovf,unf);
  291. }
  292.  
  293.  
  294. #ifdef    Faster
  295.  
  296. int_hconvolve(mask,masksize, inb, outb, r,c)
  297. int    *mask, *inb, *outb;
  298. {
  299. int    row_n,        /* row # -- position    */
  300.     v_lft_cln,    /* very left column    */
  301.     i, hrzn_p,    /* working point    */
  302.     *tmpp, *p, pc,    /* relative column pos    */
  303.     *bp=outb,
  304.     vlpix, vrpix,
  305.     sum, lweight, rweight,
  306.  
  307. last_cln= c - 1,
  308. msk2    = masksize / 2,
  309.  
  310. lftregion = msk2 > c ? c : msk2,
  311. rt_rgn_w = c - (masksize-msk2);    if (rt_rgn_w < lftregion) rt_rgn_w = 0;
  312.  
  313. for (row_n = v_lft_cln = 0; row_n < r; row_n++, v_lft_cln += c)
  314.   {    tmpp = inb + v_lft_cln;    /* current row & first column    */
  315.     vlpix = *tmpp;        /* cur row very left side pixel    */
  316.     vrpix = *(tmpp + last_cln);    /* very rlght side pix    */
  317.     for(hrzn_p = 0; hrzn_p < lftregion; hrzn_p++) {
  318.         pc = hrzn_p - msk2;
  319.         p = tmpp + pc;
  320.         /* PC moved form left half mask to right half mask    */
  321.         for(i=sum=lweight=rweight = 0; i < masksize; i++,pc++,p++)
  322.            {
  323.             if (pc < 0)
  324.                 lweight += mask[i];
  325.             else    if (pc < c)
  326.                     sum += mask[i] * *p;
  327.             else    rweight += mask[i];
  328.            }
  329.         *bp++ = sum + vlpix*lweight + vrpix*rweight;
  330.     }
  331.     for(; hrzn_p <= rt_rgn_w; hrzn_p++) {
  332.         p = tmpp + hrzn_p - msk2;
  333.         for(i = sum = 0; i < masksize;i++,p++)
  334.             sum += *p * mask[i];
  335.         *bp++ = sum;
  336.     }
  337.     for(; hrzn_p < c; hrzn_p++) {
  338.         pc = hrzn_p - msk2;
  339.         p = tmpp + pc;
  340.         for(i = sum = rweight = 0; i < masksize; i++,pc++,p++) {
  341.             if (pc > last_cln)
  342.                 rweight += mask[i];
  343.             else    sum += *p * mask[i];
  344.         }
  345.         *bp++ = sum + vrpix*rweight;
  346.     }
  347.    }
  348. }
  349.  
  350. int_vconvolve(mask, masksize, inb, outb, r, c)
  351. int    *mask, *inb, *outb;
  352. {
  353. int    row_n, hrzn_p, i,
  354.     sum, bweight, tweight,
  355.     *tmpp, *p, pr, *bp = outb,
  356.  
  357. msk2    = masksize / 2,
  358. last_row= r - 1,
  359. nd_r_p    = last_row * c,    /* last row (end row) first column    */
  360.  
  361. topregion= msk2 > r ? r : msk2,
  362. btm_rgn_h= r-(masksize-msk2); if (btm_rgn_h < topregion) btm_rgn_h = 0;
  363.  
  364. for (row_n = 0; row_n < topregion; row_n++)
  365.     for (hrzn_p = 0, tmpp = inb + (row_n-msk2)*c; hrzn_p < c; hrzn_p++)
  366.     {
  367.     p = tmpp + hrzn_p;    pr = row_n - msk2;
  368.     for(i=sum=bweight=tweight=0; i<masksize; i++,p+=c,pr++)
  369.        {
  370.         if (pr < 0)        bweight += mask[i];
  371.         else    if (pr < r)    sum += *p * mask[i];
  372.             else    tweight += mask[i];
  373.        }
  374.     *bp++ = sum + inb[hrzn_p]*bweight + inb[hrzn_p + nd_r_p]*tweight;
  375.     }
  376. for(; row_n <= btm_rgn_h; row_n++)
  377.     for(hrzn_p = 0, tmpp = inb + (row_n-msk2)*c; hrzn_p < c; hrzn_p++) {
  378.     for(i=sum= 0, p = tmpp + hrzn_p; i < masksize; i++, p+=c)
  379.         sum += *p * mask[i];
  380.     *bp++ = sum;
  381.     }
  382. for(; row_n < r; row_n++)
  383.     for(hrzn_p = 0, tmpp = inb + (row_n-msk2)*c; hrzn_p < c; hrzn_p++) {
  384.     p = tmpp + hrzn_p;
  385.     for(i=sum=tweight= 0,pr = row_n - msk2; i < masksize;i++,pr++,p+=c)
  386.         if (pr > last_row)
  387.             tweight += mask[i];
  388.         else     sum += *p * mask[i];
  389.     *bp++ = sum + tweight * inb[nd_r_p + hrzn_p];
  390.     }
  391. }
  392.  
  393. hconvolve(mask,masksize, inb, outb,r,c)
  394. Float    *mask, *inb, *outb;
  395. {
  396. int    row_n, hrzn_p, i, pc,
  397.     v_lft_cln;
  398. Float    vlpix, vrpix,
  399.     sum, lweight, rweight,
  400.     *tmpp, *p,
  401.     *bp = outb;
  402.  
  403. int
  404. last_cln    = c - 1,
  405. msk2    = masksize / 2,
  406.  
  407. lftregion= msk2 > c ? c : msk2,
  408. rt_rgn_w= c-(masksize-msk2); if (rt_rgn_w < lftregion) rt_rgn_w = 0;
  409.  
  410.   for (row_n = v_lft_cln = 0; row_n < r; row_n++, v_lft_cln += c)
  411.   {    tmpp = inb + v_lft_cln;    /* current row & first column    */
  412.     vlpix = inb[v_lft_cln];    /* cur row very left side pixel    */
  413.     vrpix = inb[v_lft_cln + c - 1];    /* very rlght side pix    */
  414.     for(hrzn_p = 0; hrzn_p < lftregion; hrzn_p++) {
  415.         pc = hrzn_p - msk2;
  416.         p = tmpp + pc;
  417.         /* PC moved form left half mask to right half mask    */
  418.         for(i=sum=lweight=rweight = 0; i < masksize; i++,pc++,p++)
  419.            {
  420.             if (pc < 0)
  421.                 lweight += mask[i];
  422.             else    if (pc < c)
  423.                     sum += mask[i] * *p;
  424.             else    rweight += mask[i];
  425.            }
  426.         *bp++ = sum + vlpix*lweight + vrpix*rweight;
  427.     }
  428.     for(; hrzn_p <= rt_rgn_w; hrzn_p++) {
  429.         p = tmpp + hrzn_p - msk2;
  430.         for(i = sum = 0; i < masksize;i++,p++)
  431.             sum += *p * mask[i];
  432.         *bp++ = sum;
  433.     }
  434.     for(; hrzn_p < c; hrzn_p++) {
  435.         pc = hrzn_p - msk2;
  436.         p = tmpp + pc;
  437.         for(i = sum = rweight = 0; i < masksize; i++,pc++,p++) {
  438.             if (pc > last_cln)
  439.                 rweight += mask[i];
  440.             else    sum += mask[i] * inb[v_lft_cln+pc];
  441.         }
  442.         *bp++ = sum + vrpix*rweight;
  443.     }
  444.    }
  445. }
  446.  
  447. vconvolve(mask,masksize,inb,outb,r,c)
  448. Float    *mask, *inb,*outb;
  449. {
  450. int    row_n, hrzn_p, i,pr;
  451. Float    sum , bweight , tweight,
  452.     *tmpp, *p,
  453.     *bp = outb;
  454.  
  455. int    last_row = r - 1,
  456. msk2    = masksize/2,
  457. nd_r_p    = last_row * c,    /* last row (end row) first column    */
  458.  
  459. topregion= msk2 > r ? r : msk2,
  460. btm_rgn_h= r-(masksize-msk2); if (btm_rgn_h < topregion) btm_rgn_h = 0;
  461.  
  462. for (row_n = 0; row_n < topregion; row_n++)
  463.     for (hrzn_p = 0, tmpp = inb + (row_n-msk2)*c; hrzn_p < c; hrzn_p++)
  464.     {
  465.     p = tmpp + hrzn_p;    pr = row_n - msk2;
  466.     for(i=sum=bweight=tweight=0; i<masksize; i++,p+=c,pr++)
  467.        {
  468.         if (pr < 0)        bweight += mask[i];
  469.         else    if (pr < r)    sum += *p * mask[i];
  470.             else    tweight += mask[i];
  471.        }
  472.     *bp++ = sum + inb[hrzn_p]*bweight + inb[hrzn_p + nd_r_p]*tweight;
  473.     }
  474. for(; row_n <= btm_rgn_h; row_n++)
  475.     for(hrzn_p = 0, tmpp = inb + (row_n-msk2)*c; hrzn_p < c; hrzn_p++) {
  476.     for(i=sum= 0, p = tmpp + hrzn_p; i < masksize; i++, p+=c)
  477.         sum += *p * mask[i];
  478.     *bp++ = sum;
  479.     }
  480. for(; row_n < r; row_n++)
  481.     for(hrzn_p = 0, tmpp = inb + (row_n-msk2)*c; hrzn_p < c; hrzn_p++) {
  482.     p = tmpp + hrzn_p;
  483.     for(i=sum=tweight= 0,pr = row_n - msk2; i < masksize;i++,pr++,p+=c)
  484.         if (pr > last_row)
  485.             tweight += mask[i];
  486.         else     sum += *p * mask[i];
  487.     *bp++ = sum + tweight * inb[nd_r_p + hrzn_p];
  488.     }
  489. }
  490.  
  491.  
  492. #else    /* easy to read and to understand    */
  493.  
  494. int_hconvolve(mask,masksize, inb, outb,r,c)
  495. int    *mask, *inb, *outb;
  496. {
  497. int    row_n,        /* row # -- position    */
  498.     v_lft_cln,    /* very left column    */
  499.     i, hrzn_p,    /* working point    */
  500.     msk2, pc,    /* relative column pos    */
  501.     *bp=outb,
  502.     vlpix, vrpix,
  503.     sum, lweight, rweight;
  504.  
  505. msk2    = masksize / 2;
  506.  
  507. for(row_n = v_lft_cln = 0; row_n < r; row_n++, v_lft_cln += c)
  508.    {
  509.     vlpix = inb[v_lft_cln];        /* very left side pixel    */
  510.     vrpix = inb[v_lft_cln + c - 1];    /* very rlght side pix    */
  511.     for(hrzn_p = 0; hrzn_p < c; hrzn_p++) {
  512.         pc = hrzn_p - msk2;
  513.     /* PC moved form left half mask to right half mask    */
  514.         for(i = sum = lweight = rweight = 0; i < masksize; i++, pc++)
  515.            {
  516.             if (pc < 0)
  517.                 lweight += mask[i];
  518.             else    if (pc < c)
  519.                     sum += mask[i] * inb[v_lft_cln + pc];
  520.             else    rweight += mask[i];
  521.            }
  522.         *bp++ = sum + vlpix*lweight + vrpix*rweight;
  523.     }
  524.     }
  525. }
  526.  
  527. int_vconvolve(mask, masksize, inb, outb, r, c)
  528. int    *mask, *inb, *outb;
  529. {
  530. int    row_n, hrzn_p, i, msk2,
  531.     last_row,
  532.     sum , lweight, tweight,
  533.     pr, tp, *bp = outb;
  534.  
  535. msk2    = masksize / 2;
  536. last_row= (r - 1) * c;
  537.  
  538. for(row_n = 0; row_n < r; row_n++)
  539.    for(hrzn_p = 0; hrzn_p < c; hrzn_p++)
  540.     {
  541.     for(i=sum=lweight=tweight=0, tp=(row_n-msk2)*c; i<masksize; i++,tp+=c)
  542.        {    pr = row_n - msk2 + i;
  543.         if (pr < 0)
  544.             lweight += mask[i];
  545.         else    if (pr < r)
  546.                 sum += mask[i] * inb[tp + hrzn_p];
  547.         else    tweight += mask[i];
  548.        }
  549.     *bp++ = sum + inb[hrzn_p]*lweight + inb[hrzn_p + last_row]*tweight;
  550.     }
  551. }
  552.  
  553. hconvolve(mask,masksize, inb, outb,r,c)
  554. Float    *mask, *inb, *outb;
  555. {
  556. int    row_n, hrzn_p, i, pc, msk2,
  557.     v_lft_cln;
  558. Float    vlpix, vrpix,
  559.     sum, lweight, rweight,
  560.     *bp = outb;
  561.  
  562. msk2    = masksize / 2;
  563.  
  564. for(row_n = v_lft_cln = 0; row_n < r; row_n++, v_lft_cln += c)
  565.    {
  566.     vlpix = inb[v_lft_cln];
  567.     vrpix = inb[v_lft_cln + c - 1];
  568.     for(hrzn_p = 0; hrzn_p < c; hrzn_p++) {
  569.         pc = hrzn_p - msk2;    /* start from over left to right */
  570.         for(i = sum = lweight = rweight = 0; i < masksize; i++, pc++) {
  571.             if (pc < 0)
  572.                 lweight += mask[i];
  573.             else    if (pc < c)
  574.                     sum += mask[i] * inb[v_lft_cln + pc];
  575.             else    rweight += mask[i];
  576.         }
  577.         *bp++ = sum + vlpix*lweight + vrpix*rweight;
  578.     }
  579.    }
  580. }
  581.  
  582. vconvolve(mask,masksize,inb,outb,r,c)
  583. Float    *mask, *inb,*outb;
  584. {
  585. int    row_n, hrzn_p,i, pr, msk2,
  586.     last_row, tp;
  587. Float    sum , bweight , tweight,
  588.     *bp = outb;
  589.  
  590. msk2    = masksize/2;
  591. last_row= (r-1)*c;
  592.  
  593. for(row_n = 0; row_n < r; row_n++)
  594.    for(hrzn_p = 0; hrzn_p < c; hrzn_p++) {
  595.     sum = bweight = tweight = 0;
  596.     for(i = 0, tp = (row_n - msk2) * c; i < masksize; i++, tp+=c) {
  597.         pr = row_n - msk2 + i;
  598.         if (pr < 0)
  599.             tweight += mask[i];
  600.         else    if (pr < r)
  601.                 sum += mask[i] * inb[tp + hrzn_p];
  602.         else    bweight += mask[i];
  603.     }
  604.     *bp++ = sum + inb[hrzn_p]*tweight + inb[hrzn_p + last_row]*bweight;
  605.    }
  606. }
  607. #endif    Faster
  608.